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STOCHASTIC NETWORKS WITH MULTIPLE STABLE POINTS 

By Nelson Antunes, Christine Fricker, Philippe Robert 
and Danielle Tibi 

Universidade do Algarve, INRIA, INRIA and Universite Paris 7 

This paper analyzes stochastic networks consisting of a set of fi- 
nite capacity sites where different classes of individuals move accord- 
ing to some routing policy. The associated Markov jump processes 
are analyzed under a thermodynamic limit regime, that is, when the 
networks have some symmetry properties and when the number of 
nodes goes to infinity. An intriguing stability property is proved: un- 
der some conditions on the parameters, it is shown that, in the limit, 
several stable equilibrium points coexist for the empirical distribu- 
tion. The key ingredient of the proof of this property is a dimension 
reduction achieved by the introduction of two energy functions and 
a convenient mapping of their local minima and saddle points. Net- 
works with a unique equilibrium point are also presented. 

1. Introduction. This paper studies the asymptotic behavior of a class 
of stochastic networks. A general description of the basic mechanisms of 
these systems is given in terms of a finite particle system or in terms of a 
queueing network. 

Description. As for some classical stochastic processes, like the zero 
range process (see Liggett [17]), one can give two alternative presentations 
for these networks. 

A particle system. It can be thought of as a set of sites where K different 
types of particles coexist. At a given site, for 1 < k < K, external particles 
of type k with mass A^ £ N arrive at rate A type k particle stays an 
exponential time with parameter 7^ at a site and then moves randomly to 
another site. A type k particle leaves the system at rate Hk- Mass constraint: 
The total mass of particles at any site must be less than C £ N, so that a 
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particle arriving at a site is accepted only if this constraint is satisfied, 
otherwise it is rejected from the system. 

A queueing network. It can be described as a set of identical finite capac- 
ity nodes where customers move from one node to another node uniformly 
chosen at random, being accepted if there is enough room and, otherwise, 
being rejected. If he is not rejected during his travel through the network, 
a customer leaves the network after his total service time. Different classes 
of customers access the network: Classes differ by their arrival rates, total 
service times, residence times at the nodes and also by the capacities they 
require at the nodes they visit. For example, a "light" customer will require 
one unit of capacity while a "heavy" customer may ask for a significant pro- 
portion of the total capacity of the node. External class k customers arrive 
at rate at any node. During their total service time, which ends at rate 
fj,l~, class k customers are transferred from a node to another one at rate 7&. 
A class k customer occupies units of capacity at any visited node, if 

this amount of capacity is not available, he is rejected. 

Large distributed networks and statistical mechanics. These stochastic 
networks have been introduced in Antunes, Fricker, Robert and Tibi [1] to 
represent the time evolution of a wireless network. Recent developments of 
mobile or sensor networks have given a strong impetus to the analysis of 
the associated mathematical models. See Borst, Hegde and Proutiere [2], 
Gupta and Kumar [11] and Kermarrec, Massoulie and Ganesh [16], for ex- 
ample. The point of view of statistical mathematics has been introduced in 
the analysis of stochastic networks some time ago by Dobrushin to study 
various aspects of queueing networks such as departure processes of queues, 
capacity regions or occupancy problems. See Karpelevich, Pechersky and 
Suhov [12] for a survey. The networks considered quite recently have a very 
large number of nodes, of the order of 10 5 -10 6 nodes in practice, this estab- 
lishes an even stronger connection with classical models of statistical physics. 
At the same time, due to the variety of possible topological structures and 
algorithms governing the behavior of these networks, new classes of math- 
ematical models are emerging. This is clearly a promising research area for 
statistical mechanics methods. 

Outline of the paper. Assuming Poisson arrivals and exponential distri- 
butions for the various service times and residence times, the time evolution 
of such a network is described by a Markov jump process with values in 
some finite (but large) state space. These associated Markov processes are, 
in general, not reversible, and little is known on the corresponding invariant 
distributions. 
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Mean-field convergence. These networks are analyzed under a thermo- 
dynamic scaling, that is, when the number of nodes N goes to infinity. It is 
shown, Section 2, that the process of the empirical distribution (Y (i)) of 
the system converges to some dynamical system (y(t)) verifying 

(1) j t y(t) = V(y(t)), t>0, 

where (V(y),y £ V(X)) is a vector field on V(X), the set of probability 
distributions on the finite set X defined by 

X = {n = (n k ) £ N K : A\n\ H V A K n K < C}. 

The dynamical system (y(t)) is therefore a limiting description of the original 
Markov process (Y N (t)). 

In general, the dimension of the state space V{X) of (y(t)) is quite large 
so that it is difficult to study this dynamical system in practice (an impor- 
tant example considered at the end of the paper is of dimension 22). The 
classification of the equilibrium points of (y(t)) with regard to the stability 
property is the main problem addressed in this paper. As it will be seen (cf. 
Proposition 5) the analysis of these points gives also insight on the limiting 
behavior of the invariant distribution of Markov process (Y N (t)). 

As a first result, it is shown that the equilibrium points are indexed by a 
"small" .ff-dimensional subset of (the number K of different classes of 
customers is usually quite small). They are identified as those elements of 
the family of probability distributions v p on X indexed by p = (p^) E . 



^ K nk 

Z \P) k=i Uk] 



where Z(p) is the partition function, for which p satisfies the fixed point 
equations 

(2) A» + TfcSn^"Wm) 

Ik + P-k 

The probability u p can also be seen as the invariant distribution of a mul- 
ticlass M/M/C/C queue. In fact the explicit limiting dynamics for the em- 
pirical distribution of the nodes is formally similar to the evolution equation 
for the probability distribution of the multiclass M/M/C/C queue with 
arrival rates Xk , service rates pk + Ik and capacity requirements , with 
the following crucial difference: the "external" arrival rates are supple- 
mented by "internal" arrival rates (corresponding to the mean arrival rates 
due to transfers from other nodes) which depend on the current state of the 
network. 

Although the equilibrium points are indexed by a subset of , a dimen- 
sion reduction of the dynamical system (y(t)) on V(X) to some dynamical 
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system of does not seem to hold. This intriguing phenomenon has also 
been noticed by Gibbens, Hunt and Kelly [10] for a different class of loss 
networks. See below. 

It is shown in Section 3 that there is a unique equilibrium point when 
all the capacity requirements of customers are equal. For arbitrary capacity 
requirements, a limiting regime of the fixed point equations with respect to 
p is also analyzed: The common capacity C of the nodes goes to infinity and 
the arrival rates are proportional to C. In this context, Theorem 2 shows 
that there is essentially one unique solution: It is shown that, if pc is a 
solution of equation (2) for capacity C then, in the limit, pc ~ i]C, where rj 
is some vector with an explicit representation in terms of the parameters of 
the network. 

A Lyapunov function. In Section 4, going back to the general case, a 
key entropy-like function g is then introduced and shown to be a Lyapunov 
function for (y(t)) so that the local minima of g on V{X) correspond to the 
stable points of (y(t)). Because of the order of magnitude of dimV(X), the 
identification of the local minima of g on V{X) is still not simple. Using this 
Lyapunov function, it is proved that, if ttn is the invariant distribution of 
the process of the empirical distributions (Y (t)), then the support of any 
limiting point of (ttn) is carried by the set of equilibrium points of (y(t)). 
In particular, when there is a unique equilibrium point y^ for (y(t)), the 
sequence of invariant distributions (ttn) converges to the Dirac mass at y^. 

A dimension reduction. A second key function <j) on the lower dimen- 
sional space is introduced in Section 5. The main result of the paper, 
Theorem 3, establishes a one to one correspondence between the local min- 
ima of g on V{X) and the local minima of <f>. The dimension reduction for 
the problem of stability of the equilibrium points is therefore achieved not 
through dynamical systems but through the energy functions g and (p. This 
result is interesting in its own right and seems to be a promising way of 
studying other classes of large networks. 

Phenomenon of bistability for (y(t)). With these results, an example of 
a network with two classes of customers and at least three zeroes for V is 
exhibited in Section 6, two of them being "stable" and the other one be- 
ing a saddle point. In this case the asymptotic dynamical system (y(t)) has 
therefore a bistability property. This suggests the following (conjectured) 
bistability property for the original process describing the state of the net- 
work: it lives for a very long time in a region corresponding to one of the 
stable points and then, due to some rare events, it reaches, via a saddle point, 
the region of another stable point and so on. This conjectured phenomenon 
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is known as metastability in statistical physics. A formal proof of this phe- 
nomenon seems to be quite difficult to obtain. The only tools available in 
this domain use either the Gibbsian characteristics of the dynamics (see 
Olivieri and Vares [19]) or at least the reversibility of the Markov process, 
Bovier [3, 4]. None of these properties holds in our case. 

Note that in Antunes et al. [1] it is proved that, for similar networks 
under Kelly's scaling, there is a unique equilibrium point. Contrary to the 
model considered here, the capacity requirement of a customer in [1] does 
not depend on his class. On the other hand, the networks analyzed here 
have a symmetrical structure: all the nodes have the same capacity and the 
routing is uniform among all the other nodes. In Antunes et al. [1], the 
routing mechanisms are quite general. 

Multiple stable equilibrium points and local dynamics in stochastic net- 
works. Asymptotic dynamical systems with multiple stable points are quite 
rare in stochastic networks. Gibbens, Hunt and Kelly [10] have shown, via 
an approximated model, that such an interesting phenomenon may occur in 
a loss network with a rerouting policy. Marbukh [18] analyzes, also through 
some approximation, the bistability properties of similar loss networks. The 
dynamics of the networks of Gibbens, Hunt and Kelly [10], Marbukh [18] 
and of this paper are local in the following sense: The interaction between 
two nodes only depends on the states of these two nodes and not on the 
state of the whole network. 

A key feature of these models is the subtle interplay between the local 
description of the dynamics and its impact on the macroscopic state of the 
network through the existence of several equilibrium points. In statistical 
physics these phenomena have been known for some time. See den Hol- 
lander [8], Olivieri and Vares [19] and the references therein for a general 
presentation of these questions. See also Bovier [3, 4] for a potential theo- 
retical approach in the case of reversible Markov processes and Catoni and 
Cerf [6] for a study of the saddle points of perturbed Markov chains. For 
the more classical setting of global dynamics, the large deviation approach 
is developed in Freidlin and Wentzell [9]. 

Phase transitions in uncontrolled loss networks. This is a related topic. It 
is known that, for some loss networks on infinite graphs, there may be several 
equilibrium distributions. If the loss network is restricted to a finite sub- 
graph T of the infinite graph Q, its equilibrium distribution 7Tjf is uniquely 
determined. It turns out that, depending on the boundary conditions on 
the sequence (ir^, T C Q), may have distinct limiting values which are 
invariant distributions for the case of the infinite graph. Significant results 
have already been obtained in this domain. See Spitzer [22], Kelly [15] for 
a survey, Zachary [24], Zachary and Ziedins [25] and Ramanan, Sengupta, 
Ziedins and Mitra [20]. 
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2. The asymptotic dynamical system. Two nodes i, j £ {1, . . . , N} of the 

network interact through the exchange of customers at rate of the order of 
1/N. Due to the symmetrical structure of the network, a stronger statement 
holds: the impact on i of all nodes different from i appears only through 
some averaged quantity. For 1 < k < K , the input rate of class k customers 
at node i from the other nodes is 



N ■ 

!<3<N,j& 

If this quantity is close to 7/ c E(A^ v fc (t)), a mean field property is said to 
hold. Note that, if the network starts from some symmetrical initial state, 
the random variables Xj^(i), j = 1, . . . , N , have the same distribution. 

Theorem 1. // 1^(0) converges weakly to z G V(X) as N tends to 
infinity, then (Y (t)) converges in the Skorohod topology to the solution 
(y(t)) of the ordinary differential equation 

(3) y'{t) = V(y{t)), 

where (y(t)) is the solution starting from y(0) = z and, for y £ V(X), the 
vector field V{y) = (V n (y),n G X) on V{X) is defined by 

K 

v n{y) = E( Afe + 7k{h,y))(yn-f k t{n k >i} - y n i{n+f k eX}) 

(4) 

+ +Vk)((n k + l)y n +f k t{ n +f k ex} ~ n k y n ), 

k=l 

where (Ik,y) = J2 m ^x m kym and f k is the kth unit vector of 



By convergence in the Skorohod topology, one means the convergence in 
distribution for Skorohod topology on the space of trajectories. 

Note that equation (4) gives the derivative dy n {t)/dt = V n {y{t)) of y n (t) 
as increasing proportionally to the difference y n -f k — yn by some factor 
Afc + 7fc(Ifc> y) which measures the speed at which nodes in state n — f k turn 
to state n (due to an arrival of some type k customer). In this factor, 7^(1^,2/) 
is added to the external arrival rate of class k customers at any node, 
and hence appears as the internal arrival rate of class k customers at any 
node. This feature characterizes the mean field property. Indeed, (Ifc, y) is the 
mean number of class k customers per node when the empirical distribution 
of the N nodes is y; so that 7fc(Ife ; y) is the mean emission rate per node of 
class k customers to the rest of the network. 
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Proof OF Theorem 1. The martingale characterization of the Markov 
jump process (Y^(t)), see Rogers and Williams [21], shows that 



M»(t)=Y n N (t)-Y»(0)- Yl n N (Y N (s),w)(w n -Y n N (s))ds 

weT(X)\{Y*(s)} 



is a martingale with respect to the natural filtration associated to the Poisson 
processes involved in the arrival processes, service times and residence times. 
By using the explicit expression of the Q-matrix Qjy, trite (and careful) 
calculations finally show that the following relation holds: 

r„ N (t) = Y„"(0) + M„ N (t) 



JU k=l \ m&X / 



(5) x (Yj?_ fk (s)t {nk > 1} - Y n N (s)t {n+fkeX} )ds 

rt K 
10 



ft K 

+ / Y(lk + »k)((n k + l)Y£ fk (s)l {n+fkeX} -n k Y^( S ))ds 

J j i 



fc=l 

t K 



+ 



J0 k=l 



From there, with a similar method as in Darling and Norris [7], it is not 
difficult to prove that if Y N (0) converges to z then: 

- the sequence (Y N (t)) of process is tight for the Skorohod topology; 

- any limit (y(t)) of (Y (t)) is continuous and satisfies the following de- 
terministic differential equation, y(0) = z and 



A 



m = e 



k=l 



Afe + 7fc E m kym.(t) 



K 



[yn-f k (t)t{n k >l} ~ yn{t)t{n+f k& X}] 



+ E(^ + Vk)[(nk + l)yn+f k (t)l{ n +f k ex} ~ n k y n (t)]. 

k=l 

This is exactly equation (4). The uniqueness of the solution of this differential 
equation implies that such a limiting point (y(t)) is necessarily unique and 
therefore that (Y N (t)) converges in distribution to (y(t)). The proposition 
is proved. □ 



The equilibrium points of the dynamical system defined by equation (4) 
are the probability distributions y £ V{X) on X such that V(y) is zero. This 
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condition can be written as follows: For n € X, 

(j2K x k + -/ k (h,y})^{n+f k ex} + (7* + Pk)n k ]j Vn 

K 

(6) = ^2(h+ik(h,y))yn-f h t {nk > 1} 

k=l 

+ (7k + fJ-k)(n k + l)y n+fk l {n+fkeX} . 

These equations are equivalent to local balance equations for the numbers 
of customers of a classical M/M/C/C queue with K classes of customers 
such that, for 1 < k < K, class k customers: 

— arrive at rate X k + 7^ (If- , y) ; 

— are served at rate 7fc + /ifc ; 

— require capacity A k . 

Consequently, (y n ) is the invariant distribution of this queue. It is well known 
(see Kelly [13], e.g.) that necessarily 

(7) yn = M^) d = Vr -,I[% n £ X, 

Z(P) ^=1 n k- 

where, for 1 < k < K, p k is the ratio of the feth arrival and service rates, 

(8) Pk 



Ik + Pk 

where (Ife,^p) is the average value of the fcth component under the proba- 
bility distribution v p on X and Z{p) is the normalization constant; in other 
words, the partition function 

K n k 

nexk=i ,Lk - 

The equilibrium points of the dynamical system (y(t)) are indexed by 
whose dimension is much smaller than V(X) thereby suggesting a possible 
simpler description of the asymptotic behavior of the network. Despite this 
quite appealing perspective, it turns out that such a dimension reduction 
cannot be achieved directly since the subset {f p ,pGlR^'} of V(X) is not left 
invariant by the dynamical system (y(t)). 

Denote by B k {p) the blocking probability of a class k customer in an 
M/M/C/C queue at equilibrium with K classes and loads pi, . . . , px , that 
is, 

1 K n h 

Z ^P> n:n+Mxh=1 n h\ 
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it is easily checked that (Ifc,f p ) = p k (l — B k (p)), so equation (8) becomes 
then 

PkPk = Afe - JkpkBk(p)- 
The following proposition summarizes these results. 



Proposition 1. The equilibrium points of the dynamical system (y(t)) 
defined by equation (3) are exactly the probability distributions v p on X , 

i Km 

W ^ = W)U^' n=MeX ' 

where Z(p) is the partition function, 

K me 

max i=i rn£ - 

and p = (pk, 1 < k < K) is a vector o/]R+ satisfying the system of equations 

(K m / K ni\ 

n:n+MXe=i n£ -' neX t=\ / 

There always exists at least one equilibrium point. 



Proof. Only the existence result has to be proved. According to equa- 
tions (7) and (8), y is an equilibrium point if and only if it is a fixed point 
of the function 

V{X) — >V(X), 

with p(y) = {pk{y)) and, for 1 < k < K, p k (y) = (\ k + -y k (h,y))/(Pk + Ik)- 
This functional being continuous on the convex compact set V(X), it nec- 
essarily has a fixed point by Brouwer's theorem. □ 



A similar situation occurs in Gibbens, Hunt and Kelly [10] where the equi- 
librium points are also indexed by the solutions p £ R+ of some fixed point 
equations and the bistability properties of the system are analyzed through 
numerical estimates. Here, a detailed study of the stability properties of the 
equilibrium points is achieved. 
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Notation. In the following, we will denote 

n K n k 

for n = (rafc) E X and p = (p^) E • The system of equations (10) can then 
be rewritten as 

3. Uniqueness results. In this section, several situations in which the 
asymptotic dynamical system (3) has a unique equilibrium point, that is, 
when the fixed point equation (10) have a unique solution, are presented. 

3.1. Networks with constant requirements. It is assumed that all classes 
of customers require the same capacity, that is, that A^ = A for k = 1, . . . , K. 
By replacing C by [C/AJ , it can be assumed that A = 1. In this case, if |n| 
denotes the sum of the coordinates of n E X , 




B\{6) can be represented as the stationary blocking probability of an M/M/C/C 
queue with one class of customers and arrival rate 9 and service rate 1. 
In this case, fixed point equation (10) are 

(11) Pk = Xk/(pk + lkB 1 (S)), k=l,...,K, 

with S = pi + • • • + pk- By summing up these equations, one gets that S is 
the solution of the equation 



H k + 7f.Bi{S) ' 

It is easily checked that S — ► B\(S) is nondecreasing and therefore that the 
above equation has a unique solution. The uniqueness of the vector (p^) 
follows from equations (11). The following proposition has been proved. 
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Proposition 2. The asymptotic dynamical system (y(t)) has a unique 
equilibrium point when capacity requirements are equal. 

In particular, when there is only one class of customers, there is a unique 
solution to equation (10). 

3.2. A limiting regime of fixed point equations. Here, the fixed point 
equations (10) are analyzed under Kelly's scaling, that is, when the capacity 
C goes to infinity and the arrival rates are proportional to C, of the order 
of X k C for the kth class. For 1 < k < K, the total service rate p k and the 
rate of residence time j k are kept fixed. It will be shown that, in this case, 
there is a unique equilibrium point. Let (M k , 1 < k < K) be a sequence of K 
independent Poisson processes with intensity 1. As usual N k {A) will denote 
the number of points of the fcth process in the subset A of R+. For C > 0, 
denote by pc = (p k (C), 1 < k < K), a solution of the fixed point equations 



\ k C = p k {C)(p k + lk J2 ^<k<K, 

\ n:n+f k £X ^ 1 n£X U ' / 



this can be rewritten as 

, 19 x x Pk(C)( , F(C7-EiIi^M([0, ft (C)])>^) ^ 

This problem is related to the limit of the loss probabilities investigated and 
solved by Kelly [14] in a general setting in terms of an optimization problem. 
The proposition below is a consequence of Kelly's result. 

Proposition 3 (Kelly's scaling). // (p k (C),l < k < K) e is such 
that 

lim Pk (C)/C = p k , l<k<K, 

C— >+oo 

with (p k ) G M.+ and 

(13) PiAi + p 2 A 2 + ---+pkAk>1, 

then, for a E N, 



c- 



nc - Sfc=i 4A4j Pk(C)]) > a) = e 
im oo P (c - Ef=i ^ fe A4([o, p fc (c)]) > o) 



-uja 



where to is the unique nonnegative solution of the equation 
(14) PiA ie ~^ + p 2 A 2 e~^ M + ■■■+ p K A K e-» A « = 1. 
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The main result of this section can now be stated. Basically, it states that, 
under Kelly's scaling, the fixed point equation (10) have a unique solution 
when the capacity gets large. 

Theorem 2. If p k > for all 1 < k < K and if for any C > the vector 
(pk(C, AC)) is any solution of equation (10) then, for 1 < k < K, 

lim Pk (C,XC)_ x k 



c^+oo C p k + 7fc - lke~ wAk 

where to > is defined as 



( ^^k + lk-lke xAk 



Proof. For 1 < k < K, the function C — > pk{C) = Pk(C, \C)/C is bounded 
by Xk/fJ-k- By taking a subsequence, it can be assumed that Pk(C) converges 
to some finite p k as C goes to infinity. Under the condition 

Al ^ + A 2 ^ + ... + A K ^>l, 

Pi P2 PK 

then necessarily p\Ai + P2A2 + • • • + PrAk > 1, otherwise one would have, 
via the law of large numbers for Poisson processes, for a > 0, 

(15) lim F(c-J2AiAfi([0,Cp i (C)])>a)=l, 

and equation (12) would then give the relation p k = \k/pk f° r 1 < k < K, 
so that 

Pi A x + p 2 A 2 + • • • + p K A K > 1 , 

contradiction. From Proposition 3 and equation (12), one gets that 

Xk = Pk(^k+lk-lke- ujAk ), l<k<K, 

were ui is the solution of equation (14) associated to {p k )- Equation (14) can 
then be rewritten as 

y X k A k e~" Ak = i 
^ Pk + lk- "i k e~" Ak 

The statement of the theorem is proved in this case. 

Now, if it is assumed that A\\\/ p\ + A2X2/P2 + ■ ■ • + ArXk / Pk < 1 then, 
since p k < Xk/ pk for all k, relation (15) holds and equation (12) finally gives 
that pk = Xk/pk, 1 < k < K . The theorem is proved. □ 
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4. An energy function on ^(X). In this section, a Lyapunov function 
is introduced. As it will be seen, it plays a key role in the analysis of the 
fixed points of the asymptotic dynamical system. Define the function g on 
the set V(X) of probability distributions on X , 

(16) g(y) = Vn log(n!j/ n ) - £ f M log Afc + lkX dx, y G V(X). 

o 

Recall that, for y G V(X), (Ik,y) = J2 m &x m kym and V{X) denotes the in- 
terior of the set V(X). 

Proposition 4. The function g is a Lyapunov function for the asymp- 
totic dynamical system (y(t)), that is, 

(V(y),Vg(y))=Y / V n (y)p-(y)<0 VyeV(X), 

^X Q y^ 

o 

and, for y E *p{X), the following assertions are equivalent: 

(a) (V(y),Vg(y))=0; 

(b) The coordinates of Vg(y) are equal; 

(c) y is an equilibrium point of (y(t)), that is, V(y) = 0. 

Proof. The vector field V(y) = (V n (y)) can be written as 

K 

v n(y) = Yl^k +7k(h, y))y n -f k t{n k >i} + {pk + ik){n k + l)y n+fk i {n+fk£X} 



k=l 



((Afc +lk(h,y))^{n+f k €X} + (Pk + lk)n k )yn] 



= Y.{F k n+h {y)-F*{y)) 
k=l 

where, for n& X, F*{y) = {^k+lk)n k yn - (\k + lk{h,y))y n -f k i{ nk >i} and 
Fn{v) = when n ^ X; note that F k = whenever n k = 0. 
For y £V(X), 

(V(y)yg(y)) = £ £ p-^n+fM ~ F niv)) 
nex k=l ° Vn 

k=ln€X \°yn-h °Vn 
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Since, for n £ X, 

dg f s - , f . v ^ , A fc +7fe(Ifc,y) 

— (y) = 1 + \og{n\y n ) - > n fc log ■ , 

dy n £rj + 7fc 

one finally gets that the relation 
(17) (V(y),Vg(y)) = £ £ log (Afc + 



fc—lneA" 



(A*fc + lk)n k y n 



holds. The quantity {V(y),'Vg(y)) is thus clearly nonpositive. On the other 
hand, for k and n such that n k > 1, 

59 / s 9g f (fJ-k + 7k)n k y n 



5j/n-/ fe V(Afe + -y k (I k , y))y n - fk 

hence (V(y),Vg(y)} is zero if and only if the coordinates of Vg(y) are equal 
and this is equivalent to the system of equations 

(Afc + Jk{h, y))Un-f k = {Hk + lk)n k y n 

for all k and n such that n& > 1, so that y is an equilibrium point of the 
asymptotic dynamical system. The equivalence (a), (b) and (c) is proved. 
□ 

Convergence of the stationary distribution. If F is some real- valued func- 
tion on M. x and y £ V(X), the functional operator associated to the Q-matrix 
Oat is given by 

tt N (F)(y)= £ n N (y,z)(F(z)-F(y)) 

zeV{X)\{y} 



£ 



£ *ky n Nl {n+fkeX} (fL + ^{e n+fk - e n )) - F(y) 
.fc=i v v J 

+ £ n k n k y n N[F[y + -^(e n _/ fc - e n )) - F(y) 



fc=i 



+ £ -^Zl n kyn(Ny m - l {m=n} ) 

\<k<K 

w / jp( „, I e "~/fe ~ e " i e m+f k — C m 

{ [ y + — n — + n 

-F(y) 
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If it is assumed that the function F is of class C 2 on R , then it is easy to 
check that the sequence (£l]y(F)(y)) converges to the following expression: 

- K K 

^2 X ky n t {n+fk( z X} (VF(y),e n+fk -e n ) + ^ fi k n k y n (VF(y), e n _ fk -e n ) 

neX U=l k=l 

+ 7fc n fc2/n?/m((VF(y),e n _/ fe -e n ) 

l<k<K 
m£X 

which is defined as Q OQ (F)(y). Moreover, by using Taylor's formula at the 
second order, this convergence is uniform with respect to y G V(X). By The- 
orem 1, one necessarily has 

n oD (F)(y) = (VF(y),V(y)), y€P(X). 

Note that this identity can also be checked directly with the above equation. 



Proposition 5. IJttn denotes the invariant probability distribution of 
(Y N (t)) on V(X), then any limiting point of (^n) is a probability distribu- 
tion carried by the equilibrium points of the asymptotic dynamical system 
(y(t)) of Theorem 1. 

In particular, if (y(t)) has a unique equilibrium point yoo, then the se- 
quence of invariant distributions (itn) converges to the Dirac mass at y c 



Joo ■ 



Proof. The set V(X) being compact, the sequence of distributions (7T/v) 
is relatively compact. Let tt be the limit of some subsequence (7Tn p )- If F is 
a function of class C 2 on M. x , then for p > 0, 

n Np (F)(y)n Np (dy) = 0. 

V(X) 

The uniform convergence of £In p (F) to (F) implies that 



0=/ n oo (F)(y)Z(dy)= (VF(y),V(y)Mdy), 
JV{X) JV(X) 

so that 7f is an invariant distribution of the (deterministic) Markov process 
associated to the infinitesimal generator ^1^. 

For t > 0, denote (temporarily) by (y(x,t)) the dynamical system starting 
from y(0) = x € V(X). Assume that there exist x € V(X) and s > such that 
y(x,s) G dV(X), that is, there exists n E X such that y n (x,s) = 0. Since 
(y n (x,t)) is nonnegative and since the function t — > y(x,t) is of class C , 
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it implies that V n (y(x,s)) = y n (x,s) = 0. This last relation, (4) defining the 
vector field (V n (y)) and the equation y n (x, s) = give that y n ± j k (x, s) = for 
any k such that n±/t£^ and consequently, by repeating the argument, all 
the coordinates of y(x, s) are null. Contradiction since y(x, s) is a probability 
distribution on X . Hence, the boundary dV(X) of V(X) cannot be reached in 
positive time by (y(x,t)). This entails, in particular, that &P(X) is negligible 
for any invariant distribution of (y(x,t)). 

For x £ V(X) and < s' < s, since the function g(y(x, •)) is of class C 1 on 
[s',s] and its derivative is (V(g)(y(x, •)), V(y(x, •))}, one has 



By integrating with respect to tt this relation, the invariance of tt for the 
process (y(x,t)) and Fubini's theorem show that 



The integrand i} 00 (g)(x) = (Vg(x),V(x)} = having a constant sign by 
Proposition 4, one deduces that 7?-almost surely, ^} 00 (g)(x) = 0. The proba- 
bility 7T is thus carried by the equilibrium points of the dynamical system. 
The proposition is proved. □ 

Asymptotic independence. In the case where (y(t)) has a unique equilib- 
rium point yoo, by using the convergence of the invariant distributions ir^ 
to the Dirac distribution 5 Voo and the fact that the coordinates of (X^ (t)) 
are exchangeable, it is easy (and quite classical) to show that for any sub- 
set / of coordinates, the random variables (X^(t),i G /) at equilibrium are 
asymptotically independent with y^ as a common limiting distribution. See 
Sznitman [23], for example. To summarize, the uniqueness of an equilibrium 
point implies that, asymptotically, the invariant distribution of the Markov 
process (X^ (t)) has a product form. 

5. A dimension reduction on M^. In this section, a function (j) on Rf is 
introduced such that p £ ffi^ is a zero of V</> if and only if the corresponding 
probability distribution v p is an equilibrium point of (y(t)). Furthermore, it 
is shown that p is a local minimum of <fi if and only if v p is a local minimum 
of g on the set of probability distributions on X . In the next section, the 



(18) 
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function <p will be used to prove the bistability of the dynamical system 
(y(t)) in some cases. 

For p = (p k ) elf, define 

K 

(19) 4>{p) = -\ogZ{p) + Y J {faPk-a k \og{ Pk )) 

k=l 

with ctfc = Afc/7fc, Pk = (ik + Pkj/lki an d Z is the partition function 

z{p)=zZ P ~ 



Proposition 6. The probability distribution v p on X is an equilibrium 
point of the asymptotic dynamical system (y(t)) if and only if\7(f)(p) = 0. 

PROOF. Remark that, for 1 < k < K, 

dZ _ p* 

so that 

hence this quantity is if and only if the fixed point equation (10) holds. 
The proposition is proved. □ 

Local minima of <j) and g. Proposition 1 has shown that an equilibrium 
point is necessarily a probability vector v p on X for some p £ . It has 
been shown that the function g defined in Section 4 decreases along any 
trajectory of the dynamical system (y(t)) by Proposition 4 so that if it 
starts in the neighborhood of a local minimum of g, ultimately it reaches 
this point. At the normal scale, that is, for a finite network, it implies that, 
with an appropriate initial state, the state of the network (X N (t)) will live 
for some (likely long) time in a subset of the states corresponding, up to a 
scaling, to this local minimum. For this reason, it is important to be able 
to distinguish stable from unstable equilibrium points of (y(t)). Due to the 
quite complicated expression defining g, it is not clear how the stability 
properties of the equilibrium points can be established directly with g. The 
function (f> plays a key role in this respect, it reduces the complexity of the 
classification of the equilibrium points according to their stability properties. 

o 

Let y £ V{X). Taylor's formula for g gives the relation, for y £ V(X), 

g(y) = g{y) + (Vg(y),y' -y)+ \y' - y)H y Jy' -y) + o(\\ y ' - y \\ 2 ), 
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where 7i y g is the Hessian matrix of g, 

and l z is the transpose of vector z. 

Propositions 1, 4 and 6 give the equivalence between: 

- y G V{X) is an equilibrium point; 
— y = v p with V<f>(p) = 0; 

- y G V(X) and (Vg(y),y> -y)=0, Vy' G V(X); 
hence the relation 

g{y>) = g(u p ) + \y> - v p )H^{y' - v p ) + o(\\y' - u p f) 

holds. 

It is assumed throughout this section that the Hessian matrix has nonzero 
eigenvalues at v p such that V0(p) = 0. Consequently, for p such that V0(p) = 

0. the probability vector v p is a local minimum of g, that is, a stable equi- 
librium point of (y(t)) if and only if the quadratic form associated to 7ig p 
satisfies the following property: 

(20) thHgPh > for all h = (h n ) with ^ h n = 0. 

It will be shown in the following theorem that relation (20) is equivalent 
to the fact that the Hessian of at p is a positive quadratic form, thereby 
establishing the dimension reduction for the problem of classification. 

Theorem 3 (Correspondence between the extrema of g and (p). 

1. A vector p G is a local minimum of the function <fi if and only if v p 
is a local minimum of the Lyapunov function g. 

2. If p is a saddle point for (ft, then u p is a saddle point for g. 

Proof. According to the above remarks, one has to study, on one hand, 
the sign of the quadratic form h — > ^Ti^h associated to g at y = u p , p G M.+ , 
on the vector space of elements h = (h n ) G IR 1 ^ such that the sum of the 
coordinates of h is 0; and on the other hand, the sign of the quadratic form 
4> at p. 

The Hessian of g and its quadratic form. It is easily checked that 

Q9 , \ -rl \- 7fc 

oy n oy m y n t^i ^k + ik\h,y) 
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The quadratic form can be expressed as 

The change of variable (h n ) — > (h n /^/y^) shows that if 

H=lh = (h n )€R x : £ Vy^ h n = 0\, 

I n£X J 

then it is enough to study the sign of the quadratic form G y on H given by 

fc=i 

where, for 1 < k < K, G i?^ is defined as 



Set 



w k d = /x i 7fc /ir =f(v^("fc - (h,y)),n e x), 

\/A fc + 7 fc (I fc ,y) 



then it is easy to check that w;^ is the orthogonal projection of vt in the 
vector space H, therefore 

K 

G y (h) = (h,h)-Y,{™i,h) 2 - 
k=i 

If W y is the subvector space of H generated by the vectors w\, 1 < k < K 
and P]v y (resp., P^rx) * s the orthogonal projection on W y (resp., on the 
orthogonal of W y ), then 

K 

G y (h) = (P w x(h),P w ±(h)) + (P Wy (h),P Wy (h)) - Y,« p w y (h)) 2 

(21) 

= (P w x(h),P w ±(h)) + G y (F Wy (h)). 

To determine the sign G y on H, it is thus enough to have the sign of G y (h) 
for h G Wy, such an element can be written as h = a\w\ + • • • + axw y K with 



( K 

G v( h ) = Yj aia i ( w h w j) ~ 

l<i,j<K V k=l 



J • 



20 N. ANTUNES, C. FRICKER, PH. ROBERT AND D. TIBI 

hence, if W y is the K x K matrix defined by W y = ((w y ,w y ), 1 < k, I < K), 

K 

(22) G y (h)= t aW y (I -W y )a ioxh = ^a k w y k . 

k=l 

The eigenvalues of the matrix W y being all real (it is symmetrical) and 
nonnegative since its associated quadratic form is nonnegative, therefore G y 
is positive on W y if and only if all the eigenvalues of W y are in the interval 
(0,1). 

The Hessian of <f> and its quadratic form. For p £ M.*f, 

d 2 <t> f , d 2 \ogZ a k 
dp k dpi dp k dpi p% 1 ' 

with (ak) = (Afe/7fe), for 1 < k,l < K. The derivatives of logZ have the 
following properties: 

a*\ d\ogZ 1 p n 



and 



#logZ 1 p» 



+ 



hence 



where 



1 V- P n 



q1 i g 2j 

-pkpi d g-j (p) = ih,v P ) ih,v P ) - (h,i,v P ) +t{k=i}(h,v P ), 



l p n 
Z(p) £ x n\ 



The quadratic form associated to at p G is given by, for a = (a k ) G M. K 
* p (a) = E ((I fc ,^)(I/,^)-(lM^p))-- + EK + (4,^))4- 

KJfc,KA- fc=l P fc 
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By using the change of variable (recall that a k = \kHk)i 



a = (a k ) -> — 

V VTfc Pk 



one gets that the sign of <j) p has the same range as the sign of *& p , where 

* p (a) = (a,a)+ £ ^ ^ 

i<k,i<K y Ajt + 7jfc(Ifc, v p ) yXi + v p ) 

x ((I k ,v p ){l h v p ) - (lk,i,v p ))a k ai 

= (a,a)- ( w k P ' w i P ) a kai, 

l<k,l<K 

with the above notation. Therefore, the sign of the quadratic form associated 
to (f> at p has the same values as the sign of ^ p (a) defined by 

(24) V p (a)= t a(I -W")a, a=(a k )eR K . 

Equations (22) and (24) show that G Up is positive on W v if and only 
if ^> p is positive on which proves assertion 1 of the theorem. Similarly, 
if p is a saddle point of <f>, equation (24) shows that the matrix W Up has 
eigenvalues in (0, 1) and in (l,+oo), so that G Up takes positive and negative 
values on W, and hence on H, v p is thus a saddle point of g. The theorem 
is proved. □ 

6. Bistability of the asymptotic dynamical system. This section gives 
an example where the asymptotic dynamical system has at least three fixed 
points: Two of them are stable and the other is a saddle point. The cor- 
responding stochastic network therefore exhibits a metastability property. 
In the limit, it suggests that its state switches from one stable point to the 
other after a long residence time. See Figure 2. The problem of estimating 
the residence time in the neighborhood of a stable point is not addressed 
here. According to examples from statistical physics, the expected value of 
this residence time should be of exponential order with respect to the size 
N of the network. For reversible Markov processes, Bovier [3, 4, 5] presents 
a potential theoretical approach to get lower and upper bounds for this 
expected value. 

A network with two classes. A simple setting is considered here: There 
are two classes of customers, K = 2, the capacity requirements are A\ = \ 
(small customers) and Ai = C (large ones) so that, at a given node, there 
may be n class 1 customers, < n < C, or only one class 2 customer. It 
is assumed that 71 = 72 = 1 and p,\ = ^2 = so that a customer leaves the 
network only when it is rejected at some node. 
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The two classes cannot coexist at a given node and, moreover, when a 
node contains class 1 customers, it has to get completely empty before ac- 
commodating a class 2 customer. Moreover, when the network is mostly filled 
with class 1 customers, the competition for capacity at each node should be 
favorable to class 1 customers, due to their large internal arrival rate (i.e., 
their arrival rate from all the other nodes). This can explain the stability of 
a state with a high density in class 1 customers. The same intuitive argu- 
ment holds for the existence of a stable state with a comparatively higher 
density in class 2 customers, though it is clear that the occurrence of this 
phenomenon should depend on the compared values of the different arrivals, 
services and transfers rates. 

Proposition 7. For a network with two classes of customers such that 
A\ = 1, A2 = C, 71 = 72 = 1, p,\ = H2 = 0, for C sufficiently large, there exist 
Ai and A2 £ M + such that the corresponding energy function <f> has at least 
one saddle point and two local minima. 

From Theorem 3, one deduces that there exists a stochastic network whose 
asymptotic dynamical system has at least two stable points. 

Proof of Proposition 7. Fix p e M.\ and choose (Ai,A 2 ) eK+ so 
that p satisfies equations (8), that is, 



fc = l,2, 



(25) A fc = p k - (I fc ,v p )=p k (l- dl Qp^ (p)j > 

by relation (23), so that v p is an equilibrium point for the limiting dynamics. 
It will be assumed for the moment that C = +00. The corresponding function 
<j) is then given by 

4>(p) = -log(p 2 + e pi ) +pi + p 2 - Ai log/31 - A 2 log,92- 
Using equation (25), one gets that 

d 2 4>, , Ai P2 e^ 



dpi Pi (P2 + ePi) 2 

P2(P2 + (1 - Pi)e pi ] 



Pi{p2 + e p1 ) 2 



and 



d 2 A2 1 

dpr p '~ pI + (p 2 + e^ >u - 

If p= {pi, pi) is chosen such that the inequality p2 < (pi — l)exp(pi) holds, 
then 



^(p)<0 and ^(p)>0. 
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Fig. 1. Function cj> with one saddle point and two stable equilibrium points. Two classes 
with Ai = 0.68, A 2 = 9.0, At = 1 and A 2 = C = 20. 




3,5e+06 4e+S6 4,5e*86 5e+B6 5,5e+B6 6e*86 6,5e*B6 7e+B6 



Fig. 2. Time evolution of the proportion of nodes without class 2 particle. Case 
N = 12000 nodes, Ai = 1, A 2 = C = 5, Ai = 0.64, A 2 = 2.71, = ^ 2 = and 71 =72 = 1- 



The constant C is now assumed to be finite and sufficiently large so that 
the above inequalities with cf> in place of <f> are satisfied, p is a saddle point 
for <f>. The function <j) is given by 

<j)(p) = - log p 2 + — T +Pl + P2 - Ailogpi - A 2 logp 2 - 
V n=0 n - / 
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The function p 2 — > 4>(p~i,P2) is convex, p 2 is a strict local minimum by 
construction and therefore a global minimum. Similarly, the function p\ — > 
cf)(pi,P2) has a strict local maximum at p±, 

M{(f)(p):p = (pi,p 2 ),Pi <Pi} < 4>{p), 

M{(j)(p):p= {pi,p 2 ),Pi <Pi} < 4>{p) =mi{(j)(p):pe A}, 

with A = {(p~i,p2) '■ P2 £ \ {0}}. Since 4>((pi,P2)) converges to +00 when 
Pi or P2 converges to or +00, one concludes that the function (ft has at least 
two local finite minima, one on each side of A. The proposition is proved. 
Figure 1 gives an example of such a situation. □ 
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